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Abstract 

We  introduce  a  generalized  Orienteering  Problem  where,  as  usual,  a  vehicle  is  routed  from  a 
prescribed  start  node,  through  a  directed  network,  to  a  prescribed  destination  node,  collecting 
rewards  at  each  node  visited,  in  order  to  maximize  the  total  reward  along  the  path.  In  our 
generalization,  transit  on  arcs  in  the  network  and  reward  collection  at  nodes  both  consume  a 
variable  amount  of  the  same  limited  resource.  We  exploit  this  resource  trade-off  through  a  spe¬ 
cialized  branch-and-bound  algorithm  that  relies  upon  partial  path  relaxation  problems  which 
often  yield  tight  bounds  and  lead  to  substantial  pruning  in  the  enumeration  tree.  We  present 
the  Smuggler  Search  Problem  as  an  important  real-world  application  of  our  generalized  Orien¬ 
teering  Problem.  Numerical  results  show  that  our  algorithm  applied  to  the  Smuggler  Search 
Problem  outperforms  standard  Mixed-Integer  Nonlinear  Programming  solvers  for  moderate  to 
large  problem  instances.  We  demonstrate  model  enhancements  that  allow  practitioners  to  rep¬ 
resent  realistic  search  planning  scenarios  by  accounting  for  multiple  heterogeneous  searchers 
and  complex  smuggler  motion. 

Keywords:  orienteering  problem;  military  operations  research;  search  and  surveillance;  route 
planning;  mixed-integer  nonlinear  programming 


1  Introduction 

We  define  the  Generalized  Orienteering  Problem  with  Resource  Dependent  Rewards  (GOP-RDR), 
which  seeks  to  route  a  vehicle  along  a  simple  path  through  a  directed  network,  between  prescribed 
start  and  end  nodes,  in  order  to  maximize  the  total  reward  along  the  path.  Rewards  are  collected  by 
the  vehicle  at  each  node  visited,  where  the  reward  level  depends  on  the  amount  of  scarce  resources 
expended.  Arcs  in  the  network  are  traversed  while  consuming  the  same  limited  resources  used  for 
reward  collection.  The  path  is  constructed  so  that  the  total  resource  expenditure  is  within  given 
limits. 

The  GOP-RDR  is  a  generalization  of  the  well-known  Orienteering  Problem  (OP)  [33].  The  OP 
and  its  multi- vehicle  extension,  the  Team  OP  (TOP),  have  wide-ranging  applicability  and  have 
been  used  to  solve  many  practical  problems  in  tourism  [24,  25,  26,  36],  sports  [3,  11,  32],  military 
operations  [16,  20],  commercial  service  and  vehicle  routing  [11,  27,  30],  and  production  [14,  18].  In 
these  problems,  node  visitation  rewards  and  arc  traversal  resource  expenditures  (arc  lengths)  are 
fixed  quantities. 
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The  GOP-RDR  generalizes  the  OP  by  allowing  node  rewards  and  arc  length  to  vary  based  on 
the  amount  of  the  resources  expended  at  each  node.  GOP-RDRs  arise  in  military,  search  and 
rescue  applications,  and  law  enforcement  operations  where  the  objective  is  to  route  a  searcher  to 
find  moving  targets  in  an  area  of  interest  (AOI)  so  that  the  total  reward  garnered  by  the  search  is 
maximized.  In  these  optimal  search  problems ,  targets  can  be  thought  of  as  carrying  some  type  of 
illicit  material.  Thus,  the  reward  garnered  by  the  search  is  related  to  the  amount  of  illicit  material 
detected.  Limited  resources  (e.g.,  time,  fuel,  etc.)  are  expended  by  the  searcher  while  performing 
search  actions  in  regions  of  interest  and  while  in  transit  between  these  regions  as  targets  move  in  the 
AOI.  We  introduce  the  Smuggler  Search  Problem  (SSP),  a  path-constrained  optimal  search  problem 
in  continuous  space  and  time  as  an  important  example  of  a  GOP-RDR.  The  SSP  deals  with  the 
high  level  decision  of  routing  search  vehicles  through  subsets  of  the  AOI  called  search  regions  in  the 
presence  of  uncertain  information  about  target  whereabouts. 

GOP-RDRs  may  also  arise  in  commercial  applications  where  a  vehicle  is  routed  to  a  number  of 
locations  in  order  to  perform  a  service.  The  distance  between  locations  can  be  represented  by  a 
travel  time,  possibly  changing  with  time-of-day  or  environmental  effects,  and  the  reward  garnered  by 
performing  the  service  at  each  location  may  be  an  increasing  function  of  time  spent  at  the  location. 
For  example,  the  Red  Cross  blood  collection  problem  [38]  is  similar  to  the  GOP-RDR  in  that  it  is 
more  beneficial  to  visit  pickup  locations  later  in  the  path  because  the  visitation  reward  increases 
with  time.  This  problem  differs  from  the  GOP-RDR  in  that,  while  rewards  at  each  node  depend 
on  time  (time  being  a  resource  consumed  in  transit  between  nodes),  the  activity  of  collecting  the 
reward  does  not  require  resource  consumption.  Moreover,  the  arc  lengths  between  nodes  do  not 
vary  in  this  problem. 

Other  generalizations  of  the  OP  have  been  considered  in  multi-objective  problems  where  rewards 
may  be  functions  of  a  number  of  attribute  scores  [24,  25,  35,  36],  and  where  the  arc  length  between 
nodes  are  determined  by  general  cost  functions  [18].  The  later  reference  introduces  a  generalized 
OP  that  arises  when  transit  resources  are  not  fixed  values  but  are  determined  by  general  resource 
expenditure  functions.  The  GOP-RDR  further  generalizes  the  OP  by  allowing  reward  collection 
to  vary  as  a  function  of  resources  expended  at  nodes.  We  are  not  aware  of  any  references  in  the 
literature  on  OPs  and  related  problems  which  consider  generalizations  of  the  node  rewards  and 
arc  lengths  at  the  same  time,  nor  have  we  encountered  any  problems  where  the  effort  of  collecting 
rewards  at  nodes  is  in  direct  competition  with  that  of  transiting  between  nodes.  The  GOP-RDR 
appears  to  be  the  first  to  consider  these  issues. 

Several  exact  algorithms  for  solving  OPs  and  TOPs  have  been  proposed  in  the  literature:  see 
[2,  4]  for  branch-and-price  references  or  [15,  19]  for  branch-and-bound  references.  Laporte  and 
Martello  [15]  describe  a  branch-and-bound  algorithm  where  fathoming  is  accomplished  by  computing 
inexpensive  upper  bounds  based  on  a  binary  knapsack  problem.  This  is  possible  because  the  arc 
lengths  and  rewards  are  fixed  values,  conditions  which  do  not  necessarily  hold  in  the  GOP-RDR. 
Lagrangian  relaxation  is  used  within  a  branch-and-bound  procedure  by  Ramesh  at  al.  [19].  They 
relax  the  budget  constraint  and  solve  the  resulting  relaxation  for  fixed  Lagrange  multipliers  using 
a  polynomial  time  degree-constrained  spanning  tree  algorithm,  a  technique  that  is  not  possible 
for  the  nonlinear  GOP-RDR.  We  present  a  branch-and-bound  algorithm  for  the  GOP-RDR  that 
capitalizes  on  the  trade-off  between  transit  and  reward  collection  resource  usage  by  solving  partial 
path  relaxation  problems  to  compute  upper  bounds.  Several  heuristics  for  solving  OPs  and  TOPs 
have  also  been  proposed  in  the  literature  [1,  6].  We  use  the  heuristic  presented  in  [6]  as  a  point  of 
departure  in  developing  a  heuristic  for  the  SSP,  where  a  simple  node  deletion  step  is  used  to  find 
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an  improving  path. 

The  remainder  of  this  article  is  organized  as  follows.  We  formulate  the  GOP-RDR  and  provide  a 
branch- and-bound  algorithm  for  obtaining  solutions  in  the  next  section.  In  Section  3  we  formulate 
the  SSP  and  describe  a  heuristic  that  is  used  to  provide  initial  solutions  to  the  branch-and-bound 
(B&B)  algorithm.  In  Section  4  we  provide  numerical  results,  comparing  branch-and-bound  solu¬ 
tions  to  solutions  obtained  by  MINLP  solvers.  Section  5  highlights  SSP  enhancements  that  allow 
practitioners  to  model  realistic  search  planning  problems.  We  conclude  with  final  remarks  in  Sec¬ 
tion  6. 


2  Formulation  and  Branch-and-Bound  Framework 

Before  formulating  the  GOP-RDR  we  begin  with  a  standard  OP  formulation  [33],  which  we  will  use 
as  a  stepping  stone  for  generalization.  Consider  a  standard  OP,  where  a  vehicle  is  routed  through 
a  transportation  network,  collecting  rewards  at  each  node.  Let  G  =  (. N ,  A)  be  the  directed  graph 
that  models  this  transportation  network,  where  N  =  {0, 1, . . . ,  n  +  1}  is  the  node  set  and  A  is  the 
arc  set.  Nodes  0  and  n  +  1  are  the  vehicle’s  home  station  and  recovery  location  respectively;  not 
necessarily  the  same  physical  location.  We  assume  that  all  arcs  incident  to  node  0  are  outbound 
arcs,  and  that  all  arcs  incident  to  node  n  +  1  are  inbound  arcs.  For  notational  convenience,  we 
define  N  =  N\{0,n  +  1}  as  the  set  of  nodes  excluding  the  home  station,  node  0,  and  the  recovery 
location,  node  n  +  1.  At  each  node  %  e  N  it  is  possible  to  collect  a  reward  qt .  Traversing  any  arc 
(i,j)  G  A  consumes  a  fixed  resource  ihJ.  Total  resource  expenditure  is  limited  by  T.  We  model  the 
vehicle  path  on  G  using  the  binary  variables  Xij,  where  xtj  takes  on  value  1  when  arc  (i,j)  is  in 
the  path,  and  0  otherwise.  We  have  the  following  OP  formulation. 

Problem  O: 


max 

X 


E& 

j&N 


S.t.  E  UjXiJ  <  T 

(ti)eA 


E  xi,j 

i:(i,j)eA 


E 

i:(j,i)eA 


X 


E  xi,j  <  !, 

i:(i,j)eA 

E  xi,j  <  l-W'l  -  i, 

(i,j)eA: 

iJeN' 

Xij  E  {0,  1}, 


3  =  0 
Vj  G  N 
j  =  n  +  1 

Vj  e  N 

\/N'  CN,N'^(b 
V(i,  j)  e  A 


(la) 

(lb) 

(lc) 

(ld) 

(le) 

(lf) 


The  objective  (la)  accumulates  rewards  along  the  path.  Constraint  (lb)  ensures  the  resources  itJ 
expended  along  the  path  do  not  exceed  the  resource  limit  T .  Constraints  (lc)  maintain  a  balanced 
network  flow  that  starts  at  the  home  station  and  ends  at  the  recovery  location.  Constraints  (Id) 
ensure  nodes  are  visited  at  most  once.  Constraints  (le)  are  the  subtour  elimination  constraints 
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proposed  in  [8],  which  are  known  to  yield  relatively  tight  linear  programming  relaxations  [29,  chap. 
!]• 

For  notational  convenience  we  will  also  use  the  auxiliary  binary  variable  i/j,  which  is  uniquely 
determined  by  variables  xhj.  Variables  yj  take  on  value  1  when  node  j  is  in  the  path,  and  0 
otherwise;  i.e., 

l/o  =  l,  Vj=  xid,\/jeN\{  0}  (2) 

We  denote  by  x  the  vector  of  path  variables  {xt)J  :  ( i,j )  G  A}.  We  denote  by  y  the  vector  of  node 
visitation  variables  {yj  :  j  G  N}  and  represent  (2)  with  the  expression  y  =  Tx  for  an  appropriately 
selected  matrix  T.  We  define  X  as  the  set  of  paths  that  satisfy  (lc),  (Id),  (le),  and  (If). 

We  relax  O  to  construct  a  GOP-RDR  as  follows.  A  visit  to  any  node  i  G  {1, . . . ,  n}  is  rewarded 
at  the  expense  of  consuming  r  dwell  resources  d,  G  3ftr.  Similarly,  transit  resources  tl:J  G  3ftr  are 
consumed  when  traveling  directly  from  node  i  to  node  j.  Resources  may  represent,  for  example, 
various  consumables  such  as  time,  fuel,  and/or  money.  For  notational  convenience,  we  also  include 
the  auxiliary  resource  status  variable  a,  G  3ftr.  This  variable  is  used  to  track  the  accumulation  of 
resources  expended  along  the  path.  Let  a,  d  G  sRr(n+2);  and  t  G  -RdA  denote  vectors  of  resource 
variables;  for  example  a  =  (aj,  a[  , . . . ,  a^+1)T.  A  vehicle  may  collect  rewards  according  the  utility 
function  /(d )  :  3fJr(n+2)  h »  3ft.  We  assume  without  loss  of  generality  that  no  reward  is  possible  at 
nodes  0  and  n+ 1.  We  assume  that  /  is  a  concave  utility  function,  where  /( 0)  =  0.  The  vehicle  path 
through  G  must  obey  rj  resource  expenditure  laws  on  each  arc  (i,j)  G  A  denoted  by  the  functions 
hij(ai,  a,j ,  di ,  dj,tij)  :  3 ft5r  (->•  3ft^.  The  resource  expenditure  laws  account  for  applications  where  arc 
lengths  are  allowed  to  vary  as  resources.  In  the  static  network  considered  in  O,  ht,3  =  \tltJ  —tij ||  for 
fixed  arc  lengths  tjj,  but  in  Section  3  arc  lengths  are  not  fixed  because  the  nodes  of  the  network  are 
in  motion.  The  vehicle  path  must  also  obey  resource  expenditure  laws  at  each  node  i  G  N  denoted 
by  the  functions  gi(at,  d,)  :  3ft2r  (->•  3ft7  and  m^a^d,;)  :  3ft2r  i-»  3ftM.  We  assume  that  functions  hhJ 
and  (jj  are  convex,  and  functions  m,j  are  affine.  The  vehicle  path  must  be  such  that  total  resource 
expenditure  stays  within  the  resource  limits  defined  by  T  G  3ftr.  Let  the  matrix  Y  G  sftr(n+2)xr(r!+2) 
be  the  diagonal  matrix  diag(7/oO,  y{\_,  y2l, . . .  ,ynl,yn+iO),  where  0, 1  G  3ftr  are  a  vectors  of  0s  and 
Is  respectively.  We  note  that  the  expression  Y d  simply  returns  the  dwell  resource  vector  associated 
with  reward  collection  nodes  in  the  path  x.  We  now  state  the  GOP-RDR. 

Problem  P: 

Sets 

N 
A 
X 

Parameters 

T  G  3ftr  resource  expenditure  limits 

r  G  3ft(«+2)xIA  path-to-node  visitation  mapping  matrix  representing  (2) 


nodes:  i,j  G  {0, 1, . . . ,  n  +  1} 
arcs 

paths  that  satisfy  (lc),  (Id),  (le),  and  (If) 
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Functions 

/ :  sr(n+2)  >->  3? 

hij  :  9tf5r  ^  W 
gt'j  :  9?2r  ^  3T 
nnj  :  5ft2r  ^  W 


concave  reward  collection  objective  function 
convex  resource  expenditure  law  functions 
convex  node  resource  expenditure  law  functions 
affine  node  resource  expenditure  law  functions 


Variables 

a, 


di 


tij 


X 


h3 


Vi 


node  i  auxiliary  resource  variable 
node  i  dwell  resource  variable 
arc  (i,j)  transit  resource  variable 
arc  (i,j)  binary  path  variable 
node  i  binary  visitation  variable; 

Y  is  a  diagonal  matrix  of  these  variables 


Formulation 


max  f(Yd) 

a,d,t,x,y 

s.t.  dj ,  dj ,  t^Xij  ^  0, 

S  dj  +  tij  <  T 

j£N  ( i,j)eA 

9j{aj,  dj)  <  0, 

'".iioj-  dj)  =  0, 

a,j,  dj  >  0, 

tij  >  0, 

y  =  Yx 

x  G  X 


(3a) 

V(i,j)  e  A 

(3b) 

(3c) 

V]  e  N 

(3d) 

Vj  e  N 

(3e) 

Vj  e  N 

(3f) 

V(i,  j)  e  A 

(3g) 

(3h) 

(3i) 

The  objective  (3a)  maximizes  the  reward  collected  along  the  path.  Constraints  (3b)  enforce  resource 
expenditure  laws  on  each  arc.  Constraint  (3c)  ensures  that  total  resource  expenditure  is  within 
the  prescribed  limits.  Constraints  (3d)  and  (3e)  enforce  resource  expenditure  laws  at  each  node. 
Concavity  of  /  makes  it  desirable  consume  resources  d.  Constraints  (3c)  make  it  undesirable 
consume  resources  t.  However,  d  and  t ,  along  with  a,  are  related  through  constraints  (3b)  so  it 
may  not  be  possible  to  consume  t  =  0  transit  resources.  Observe  that  when  node  j  is  not  in  the 
path  x,  (3d)  and  (3e)  are  vacuous  because  aj  and  dj  can  be  chosen  arbitrarily  to  satisfy  these 
constraints  provided  total  resource  expenditure  (3c)  is  not  exceeded.  Constraint  (3c)  makes  dj  =  0 
desirable  in  this  situation  because  a  higher  reward  is  obtained  by  consuming  more  dwell  resources 
at  visited  nodes.  We  assume  that  dj  =  0  is  always  feasible.  Similarly,  when  arc  (i,j)  is  not  in  the 
path  x ,  constraint  (3b)  is  inactive  and  the  resource  constraint  (3c)  forces  ttJ  =  0.  Constraints  (3f) 
and  (3g)  require  nonnegative  resource  expenditure.  Binary  visitation  and  path  variables  are  set  by 
(3h)  and  (3i)  respectively. 

This  MINLP  as  posed  has  a  non-convex  continuous  relaxation.  When  this  is  the  case  most 
MINLP  solvers  such  as  DICOPT  [12]  or  BONMIN  [34]  provide  no  guarantees  of  finding  globally 
optimal  solutions.  The  problem  can  be  convexihed  with  a  Big-M  reformulation.  For  example,  (3b) 
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can  be  reformulated  as 


li;_i(ai.ardi.dj.tiJ)  <M(  1  -  xitj),  V(i,  j)  G  A, 

for  M  sufficiently  large.  We  discuss  Big-M  reformulation  and  show  numerical  results  in  the  context 
of  the  SSP  in  the  sequel.  Another  approach  is  to  use  a  B&B-based  MINLP  solver  such  as  BARON 
that  uses  convexifying  techniques  at  each  node  of  a  B&B  enumeration  tree  to  obtain  globally  optimal 
solutions  [22,  28].  This  approach  is  not  pursued  here,  rather  we  use  the  underlying  structure  of  P 
as  a  basis  for  computing  solutions. 

With  these  matters  in  mind,  we  now  proceed  to  describe  a  B&B  approach  that  utilizes  convex 
relaxation  problems,  avoids  Big-M  reformulations,  and  capitalizes  on  the  underlying  structure  of 
P  as  a  basis  for  branching  and  pruning.  We  introduce  the  notation  G(x )  =  (N(x),  A(x)),  where 
N(x )  =  {j  G  N  :  ijj  =  1; y0,  =  1  ,y  =  Tx}  and  A{x)  =  {(i,  j)  G  A  :  xl)3  =  1}.  For  any  path  x  G  X, 
P  can  be  expressed  as  the  following  convex  NLP. 

Problem  P(jc): 


max 

a,d,t 

f(Yd) 

(4a) 

S.t. 

CD 

VI 

TS 

<3 

■J2 

V(i,  j)  e  A(x) 

(4b) 

S  dj  +  JO  tij  <  T 

j€N{x)  (i,j)eA(x) 

(4c) 

9j(aj,  dj)  <  0, 

Vj  G  N(x) 

(4d) 

o' 

■V~J 

Vj  G  N(x) 

(4e) 

o' 

Al 

e 

Vj  G  N(x) 

(4f) 

tij  >  0, 

V(i,  j)  e  A(x) 

(4g) 

When  x  is  fixed,  y  can  be  computed  by  (2).  Variables  a3  and  dj  corresponding  to  unvisited  nodes 
are  removed  from  the  problem.  Similarly,  variables  th3  corresponding  to  arcs  not  traversed  are 
eliminated.  The  resulting  convex  NLP  in  the  remaining  variables  a ,  d,  and  t  is  efficiently  solved  by 
standard  NLP  solvers  such  as  MINOS  [17].  If  we  are  able  to  enumerate  all  possible  paths  x  G  X 
and  solve  the  associated  fixed  path  NLP,  we  are  assured  to  find  the  optimal  solution  to  P. 

Observing  that  since  d  contributes  to  reward  collection  and  t  consumes  resources  without  reward, 
nonzero  values  of  d  and  small  values  of  t  are  always  desired.  We  construct  the  matrix  I  by  taking 
an  r(n  +  2)  x  r(n  +  2)  identity  matrix,  and  setting  the  first  r  diagonal  entries  and  the  last  r  diagonal 
entries  to  zero.  The  expression  Id  returns  the  vector  of  dwell  resources,  setting  home  and  recovery 
dwell  resources  to  zero.  We  define  StJ  G  as  the  smallest  possible  resource  expenditure  between 
node  i  and  node  j.  If  we  consider  P(cc)  and  allow  reward  to  be  collected  at  every  node  with  no 
transit  resource  expenditure  on  any  arc,  we  obtain  the  following  relaxed  NLP. 

Problem  RP(O): 


max 

a,d 

f(Id ) 

(5a) 

S.t. 

Y.dj<T- 

j£N 

1 

O 

"3  oiS 

min  (<5j  „+i} 

j£N: 

[j,n+l)eA 

(5b) 

(3d),  (3e),  and  (3f) 
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Observe  that  the  resource  limit  decrement  on  the  right  hand  side  of  (5b) 

mm  {Soj}  +  min  {^,n+i}  (6) 

(OJ)eA  {j,n+\)eA 

is  a  lower  bound  on  t0j  +  tj,n+i,Vj  G  N.  A  path  that  visits  any  nonempty  subset  of  nodes 
N'  C  N  :  N'  0  consumes  at  least  (6)  transit  resources.  RP(O)  is  clearly  a  relaxation  of  P. 
P  is  obtained  by  adding  constraints  (3b),  (3g),  (3h),  and  (3i)  to  RP(O),  while  restricting  reward 
collection  (3a)  to  nodes  in  the  path,  and  incurring  a  transit  resource  expenditure  in  (3c)  that  is  no 
less  than  (6).  Denoting  the  optimal  objective  function  values  of  P.  P(cc),  and  RP(O)  by  Z*  1  Z(x)*, 
and  Z(0)*,  respectively,  we  state  the  following  result. 

Proposition  1.  Z(0)*  >  Z*  >  Z(x)*,\/x  G  X 

Proof.  The  result  follows  from  the  fact  that  RP(O)  is  a  relaxation  of  P  and  that  Z*  =  max  Z(x)*. 

□ 

In  order  to  obtain  useful  bounds  on  P.  we  introduce  the  notion  of  a  partial  path.  We  define  a 
partial  path  xg  to  be  the  binary  vector  satisfying  constraints  (Id),  (le),  and  (If),  while  constraints 
(lc)  are  satisfied  for  all  nodes  except  the  recovery  location  n  +  1  and  the  last  node  £  visited.  We 
define  the  indicator  parameter  I g  that  takes  on  value  1  when  l  —  n  +  1,  and  0  otherwise. 

For  any  partial  path  xg  we  have  the  following  convex  partial  path  relaxation  NLP. 

Problem  RP(3b): 


max 

a.d.t 

f(Id) 

(7a) 

S.t. 

kg  ]  ( (ij .  dj ,  dj ,  dj ,  t  j.,j )  L  0, 

V(z,  j)  G  A(xg) 

(7b) 

dj  +  tij  <  T  —  (1  —  Ie)6e,n+i 

jeN  ( i,j)eA(£e ) 

(7c) 

o' 

tT 

Vj  G  N\N(xg) 

(7d) 

tij  >  o, 

V(z,  j)  G  A(xg) 

(7c) 

(3d),  (3e),  and  (3f) 

We  note  that  when  Ig  =  1,  the  path  is  complete,  and  P(:r)  and  RP(ib)  are  equivalent.  Conversely, 
when  Ig  =  0,  the  path  is  a  partial  path,  (7d)  is  inactive  and  resources  associated  with  unvisited 
nodes  are  allowed  to  take  on  nonzero  values,  and  the  right  hand  side  of  (7c)  is  decremented  by 

8g,n+ 1- 

Consider  extending  the  partial  path  Xg  by  adding  any  arc  (£,  k)  G  A  to  the  path  as  shown  in 
Figure  1.  In  this  depiction  the  minimum  resource  expenditure  <5^n+i  in  the  partial  path  Xg  is  no 
larger  than  the  transit  resource  expenditure  tg ^  +  8g^n+\  in  the  partial  path  x^.  Next  we  show  a 
result  that  supports  building  successive  restrictions  of  RP(ab)  by  adding  to  the  partial  path.  Let 
Z(xg)*  be  the  optimal  objective  function  value  of  RP(ab). 

Proposition  2.  Z{x.g)*  >  Z(x^)*,\/k  :  (£,k)  G  A 

Proof.  Observe  that  if  t  —  n  +  1,  then  {k  :  (£,k)  G  A}  =  0.  Suppose  i  G  N\{n  +  1}.  Adding 
node  k  and  arc  (£,  k)  to  the  partial  path  adds  a  block  of  constraints  to  (7b)  and  variable  tg)k  in  (7c) 
and  (7e).  Since  the  increase  in  resource  expenditure  along  the  new  partial  path  is  at  least  8g,n+ 1, 
RP(^jt)  is  a  restriction  of  RP(ib)  and  the  result  follows.  □ 
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Figure  1:  Partial  path  Xe  (left).  Partial  path  Xf.  (right).  Transit  resource  expenditure  along  partial 
path  Xi  is  no  larger  than  transit  resource  expenditure  along  path  x^. 

We  present  a  GOP-RDR  B&B  framework  that  begins  at  the  home  station  and  forms  partial 
paths  by  adding  nodes  to  a  path  sequentially,  solving  restrictions  of  R,P(ap)  along  the  way.  We 
denote  the  set  X  as  the  set  of  all  possible  partial  paths  x .  Let  l  G  {0, 1, . . . ,  n  +  1}  denote  the  level 
of  the  B&B  enumeration  tree,  and  let  &  C  N  be  the  set  of  nodes  yet  to  be  considered  at  level  l. 
We  define  the  set  C  X  x  N  to  be  a  subset  of  partial  path  and  B&B  enumeration  tree  level  pairs. 
Let  e  >  0  be  the  absolute  optimality  gap  stopping  tolerance. 

Algorithm  B&B: 

1.  Initialize  i  =  0;  x*  =  xe  =  0;  lower  bound  LB  —  0;  l  —  0;  L0  —  0;  —  N,  \/k  —  1, . . . ,  n  +  1; 

and  =  {(0,0)}.  Solve  RP(0).  If  RP(0)  is  infeasible,  then  stop;  P  is  infeasible.  Otherwise, 
initialize  upper  bound  UB  =  Z(0)*. 

2.  If  UB  —  LB  <  e,  then  stop  and  return  x* .  Otherwise,  choose  (xe,  l )  G  0.  Add  node  j  to  partial 
path  xe  to  form  the  extended  partial  path  Xj  G  X  that  contains  arc  (£,  j).  Add  {( Xj ,  l  +  1)} 
to  O.  Remove  j  from  &+i.  Solve  RP(*j). 

3.  If  j  —  n  +  1  and  LB  <  Z(xj)*  and  RP(ij)  is  feasible,  then  set  LB  =  Z(ibj)*  and  x*  =  Xj  as 
the  best  complete  path  found  thus  far. 

4.  If  Ll+i  =  0  then  set  UB  =  rna x{Z*+1,  Zf+1},  where  Z*+1  is  the  largest  value  Z(x)*  of  all 
complete  paths  explored  thusfar  and  Zl*+1  is  the  largest  value  Z(£j)*  on  level  l  + 1  of  the  B&B 
enumeration  tree. 

5.  If  Z(ay*  <  LB  or  if  RP(*j)  is  infeasible,  then  fathom  partial  path  x3  by  removing  from  O 
all  elements  (x-j,l),  where  x3  is  a  subpath  of  x  -  and  l  >  l. 

6.  Return  to  step  2. 

This  algorithm  can  be  accelerated  by  obtaining  an  initial  feasible  solution  x  that  produces  a 
better  lower  bound  in  step  1,  thereby  allowing  fathoming  in  step  5  to  occur  more  rapidly.  To  this 
end,  we  provide  a  specialized  heuristic  for  the  SSP  in  Section  3.3.  We  do  not  prescribe  the  nature  of 
branching  to  be  performed  step  2.  Numerical  results  discussed  in  Section  4  use  depth-first-search, 
but  other  branching  strategies  can  also  be  used  (see,  for  example  [19]  or  [29,  ch.  2]). 
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3  The  Smuggler  Search  Problem 


The  Smuggler  Search  Problem  (SSP)  is  special  case  of  a  GOP-RDR  that  arises  in  challenging  real- 
world  search  operations.  This  work  is  motivated  by  ongoing  efforts  to  detect  and  interdict  the 
flow  of  illicit  traffickers  in  international  waters.  To  accomplish  this  mission  coalition  forces  strive  to 
employ  a  limited  number  of  search  assets  as  effectively  as  possible,  under  strict  resource  constraints, 
as  they  respond  to  uncertain  estimates  of  how  illicit  traffickers  move  in  the  AOl. 

OPs  and  similar  models  have  been  used  to  solve  path-constrained  optimal  search  problems. 
Path-constrained  optimal  search  problems  are  known  to  be  NP-hard  [31].  Many  approaches  which 
focus  on  discrete  space  and  time  models  can  be  found  in  the  literature  [7,  10,  13,  21,  23].  We 
proceed  to  formulate  the  SSP,  a  novel  path-constrained  optimal  search  model  in  continuous  space 
and  time. 

Consider  a  planning  scenario  where  a  search  vehicle  is  to  be  routed  throughout  an  AOI  to 
detect  multiple  moving  targets  and  within  a  D-hour  long  mission  execution  period.  We  assume 
that  the  targets  are  in  linear  motion  with  a  constant  speed,  independent  of  the  search  effort.  The 
search  controller  has,  based  on  planning  factors  and  intelligence  estimates,  the  information  listed 
in  Table  1. 


Maximum  cruise  speed  of  the  searcher  while  in  transit 

Speed  of  the  searcher  while  performing  search  actions 

Searcher  sensor  sweep  width 

Searcher  endurance  time  limit 

Scenario  time  limit 

Number  of  targets 

Speed  of  target  j 


V 

V 
W 
T 
D 
n 
Uj 


Table  1:  Known  data  available  to  the  search  controller 


Suppose  that  there  is  uncertainty  with  regard  to  where  and  when  each  target  departs,  as  well 
as  the  value  of  detecting  the  target,  but  the  nature  of  the  intelligence  allows  the  search  controller 
to  estimate  these  values  within  some  range  of  uncertainty.  These  data  are  listed  in  Table  2. 


Expected  departure  time  of  target  j 
Time  uncertainty  range  of  target  j 
Expected  departure  location  of  target  j 
Expected  arrival  location  of  target  j 
Departure/arrival  location  uncertainty  range  of  target  j 
Expected  value  of  detecting  target  j 


Pj 

Pj 

Pj 

Qj 


Table  2:  Uncertainty  data  derived  by  the  search  controller 


Based  on  these  data,  the  latest  departure  time  for  each  target  j  can  be  calculated  as 


^min 

tj 


U  +  \fy 


Similarly,  the  earliest  arrival  time  for  each  target  j  can  be  calculated  as 


.max  _ 

j  ~ 


2  Tj+Ujllf>j  Pj 
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Lastly,  we  can  calculate  the  velocity  vector  of  target  j,  Uj  as  a  function  of  speed,  and  expected 
departure  and  arrival  locations. 

We  model  search  within  each  target’s  region  of  uncertainty  (search  region)  using  a  random 
search  law  with  known  sensor  sweep  width  W;  for  details  on  random  search  models  see  [37,  ch. 
2],  We  assume  that  search  actions  cannot  be  conducted  for  more  than  one  target  at  the  same 
time.  This  reflects  the  operational  setting  where  the  searcher  is  seeking  out  a  specific  target  looking 
for  characteristics  outlined  in  intelligence  reports.  Thus,  in  the  event  that  search  regions  overlap 
in  space  and  time,  the  searcher  cannot  receive  additional  reward  for  searching  for  more  than  one 
target  at  a  time.  Given  all  the  available  information  the  search  controller  wishes  to  route  the  search 
vehicle  through  the  AOI  in  order  to  maximize  the  expected  value  of  the  search  effort.  We  model 
this  as  routing  a  vehicle  across  a  transportation  network  G  =  ( N,A ),  where  nodes  are  defined  by 
the  search  regions  and  arcs  are  defined  by  the  searchers  transit  between  each  pair  of  moving  search 
regions. 


3.1  Formulation 

Since  targets  are  in  linear  motion,  the  distance  required  to  travel  directly  between  each  (i,j)  pair 
of  search  regions  can  be  computed  as  a  function  of  time  using  Euclidean  distance  calculations.  We 
proceed  under  the  assumption  that  the  path  of  the  searcher  is  through  the  center  of  each  search 
region.  We  also  assume  that  the  error  between  Euclidean  distance  and  great  circle  distance  is 
small  relative  to  the  size  of  the  search  regions.  Suppose  a  searcher  is  searching  for  target  in  some 
predefined  order  and  that  the  searcher  has  just  completed  searching  region  i.  If  a*  represents  the 
time  the  searcher  began  searching  region  i  and  di  represents  the  duration  of  the  search  in  region  i, 
then  we  can  compute  the  current  position  p{  of  the  searcher  as  pt  =  p*  +  (cq  +  di  —  Ti)Ui.  We  assume 
that  di  >  rt ,  Vi ;  the  searcher  will  never  arrive  to  search  a  target  that  has  not  departed.  Suppose  the 
searcher  is  next  routed  to  region  j,  and  that  the  transit  time  from  region  i  to  region  j  is  denoted  as 
tij.  The  position  of  region  j  at  the  moment  the  searcher  arrives  is  pj  =  pj  +  (a*  +  d*  +  t,ltJ  —  Tj)uj. 
We  can  now  relate  the  distance  between  region  i  and  region  j  to  the  distance  the  searcher  can  travel 
in  the  same  amount  of  time  | \p,.  —  Pj\  j  =  Vtij.  We  can  relax  this  relationship  by  recognizing  that 
the  searcher  does  not  always  have  to  travel  at  maximum  cruise  speed.  The  searcher  could  choose  to 
travel  slower.  Thus,  \pt  —  Pj  ||  <  Wjj,  which  is  second-order  cone  constraint  in  the  time  resources 
a,  d ,  and  t  of  the  form  (3b). 

Drawing  from  search  theory  [37,  ch.  2]  we  define  the  detection  rate  in  search  region  j,  aj,  as 


=  _WV^ 

'  TjPjUj' 


(8) 


We  assume  that  the  searcher  speed  V  is  much  greater  than  the  target  speeds  Uj,\/j  E  N.  From  the 
searcher’s  perspective,  within  each  search  region,  the  target  is  essentially  stationary.  It  is  possible  to 
model  the  problem  where  this  does  not  hold  [37,  sec.  6-1].  However,  this  assumption  approximately 
holds  in  our  SSP  model,  where  we  consider  search  aircraft  and  surface  (e.g.,  boats)  smugglers. 

This  problem  can  be  formulated  as  the  following  MINLP,  which  is  a  special  case  of  P. 
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Problem  SSP: 


max 

a,d,t,x,y 

S.t. 


E  Qj  (1  -  exp  {-ajdjyj})  (9a) 

j&N 

(|| Pi  +  (a*  +  di  -  Ti)ui  -  Pj  -  (a,i  +  di  +  titj  -  Tj)Uj\\ 

...  -  Vtij)xij  <  0,  V(i,  j)  G  A  (9b) 

(a*  +  di  +  tid  -  dj)  xi>:j  <  0,  V(i,  j)  G  A  (9c) 

Edj+  E  (9d) 

jeiv  (i,j)eA 

E  dj  +  E  b.j  <  -D  (9e) 

jew  (i,j)eA 

dj  >  r™in,  Vj  G  iV  (9f) 

dj  +  dj  <  rj7iax,  Vj  G  iV  (9g) 

do  =  0  (9h) 

dn+1  =  o  (9i) 

o>ji  dj  >  0,  Vj  <E  N  (9j) 

Uj>  0,  V(i,j)Gd  (9k) 

y  =  Tx  (91) 

x  G  X  (9m) 


The  objective  (9a)  is  to  maximize  the  expected  value  of  the  search  effort.  (9b)  ensures  that  the 
distance  between  search  region  i  and  search  region  j  obeys  Pythagorean’s  Theorem.  Constraints 
(9c)  propagate  arrival  times  a  forward  in  time  as  arcs  are  traversed.  (9b)  and  (9c)  correspond  to 
(3b)  in  P.  Constraints  (9d)  and  (9e),  corresponding  to  (3c)  in  P.  ensure  that  the  plan  does  not 
exceed  resource  limits  T  and  D  respectively.  Note  that  the  left  summation  in  (9d)  is  over  the  set  of 
nodes  not  including  the  home  station  and  recovery  location  N.  This  may  appear  to  be  inconsistent 
with  (3c),  however  in  SSP  we  could  equivalently  model  two  dj  terms  for  each  node.  One  retains  the 
correct  dwell  resources  at  all  nodes,  and  the  other  is  nearly  a  copy  but  consumes  zero  dwell  resources 
at  nodes  0  and  n+  1.  We  choose  the  more  compact  formulation  here.  Constraints  (9f)  require  that 
the  vehicle  be  routed  to  search  regions  only  after  the  target  has  departed.  Similarly,  constraints 
(9g)  preclude  searching  in  a  region  after  the  target  has  arrived.  (9h)  and  (9i)  are  included  for 
completeness.  They  require  respectively  that  time  resources  start  at  0  and  that  the  mission  ends 
upon  recovery.  (9f),  (9g),  (9h)  and  (9i)  correspond  to  (3d)  in  P.  Constraints  (9j),  (9k),  (91),  and 
(9m)  are  as  discussed  previously. 

The  dwell-to-transit  resource  trade-off  that  underlies  SSP  can  be  observed  as  follows.  Consider 
the  objective  (9a),  and  constraints  (9b)  and  (9d).  Since  the  objective  (9a)  is  monotonically  increas¬ 
ing  in  d ,  large  dwell  time  is  desirable.  At  the  same  time,  since  constraint  (9d)  limits  the  searcher’s 
flying  time,  small  values  of  t  are  desirable.  However,  due  to  constraints  (9b),  going  from  one  search 
region  to  the  next  consumes  transit  time.  This  is  the  trade-off.  Dwell  time  in  a  search  region  cannot 
be  consumed  without  also  consuming  transit  time  in  order  to  get  to  said  search  region. 

We  define  an  NLP  analogous  to  P(a;)  which,  for  any  path  x  G  X,  provides  the  optimal  time 
resource  expenditure.  We  arrive  at  this  problem  by  fixing  x  and  y,  and  retaining  from  SSP  only 
the  interesting  constraints  and  objective  function  terms.  Recall  that  N(x)  =  {j  G  N  :  yj  =  l]yo,= 
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1,2/  =  Tcc}  and  A(x)  =  {( i,j )  G  A  :  x%t]  =  1}.  Additionally,  we  define  the  set  of  search  regions  in 
the  path  N(x )  =  N(x)\{0,n  +  1}. 

Problem  SSP  (a;): 

max  E  Qj  (1  —  exP  {~ajdj}) 

a  .(i  t  .  . 

jeN(x) 

s.t.  1 1  Pi  +  (a*  +  di  -  Ti)ui  -  Pj  -  ( a.i  +  di  +  titj  -  Tj)Uj\\  -  Vtid  <  0,  V(i,  j)  G  A(x) 

cti  +  di  +  tij  —  cij  <  0,  V(i,  j)  G  A(x) 

Y  dj  +  E  ti,j  —  T 

jeN(x)  *) 


E  dj  +  Y  ti,j  -  D 

j€N{x)  ( i,j)eA(x ) 


n  .  >  rr-min 

Uj  ^  , 

Vj  G  N(x) 

dj  +  dj  <  T™aX, 

Vj  G  N(x) 

dj ,  dj  >  0, 

Vj  G  N(x) 

U,j  >  0, 

V(i,j)  G  A(x) 

(9h)  and  (9i) 

After  a  Big-M  reformulation,  SSP  can  be  equivalently  stated  as 

an  MINLP  with  a  convex 

continuous  relaxation.  Let  if/A  be  a  number  that  is  always  greater  than  the  distance  between 
region  i  and  region  j.  Let  be  a  number  that  is  always  greater  than  the  time  required  for  the 
searcher  to  travel  between  region  i  and  region  j.  Let  M?  be  a  number  that  is  always  greater  than 


the  search  time  in  region  j. 

Problem  SSPM 

max 

a,d,t,x,y 

E  Qj  (1  -  exp  {-ajdj}) 

j&N 

(1  la.) 

S.t.  |  Pj  + 

(ttj  T  dj  Tj  j U -j  Pj  (di  -j-  dj  -j-  f  jj  Tj  j Uj  | 

...  -  Vtij  <  Mgj(  1  -  Xij), 

V(i,  j)  e  A 

(lib) 

di  -j-  di  -|-  ti  j  dj  A  Afjj  (1  (T  jy )  > 

V(i,j)  G  A 

(11c) 

dj  <  MPyj, 

Vj  G  N 

(lid) 

Ydj  +  Y  t 

jeN  (ti)eA 

(lie) 

Y  dj  +  Y  ti,j  <  D 

j£N  ( i,j)eA 

(Ilf) 

(9f)  through  (9m) 

The  nonlinear  interactions  between  the  binary  variables  and  the  continuous  variables  in  SSP 
are  modeled  with  a  Big-M  on  the  right  hand  sides  of  (lib),  (11c)  and  (lid).  Constraint  (lid) 
requires  that  search  duration  be  zero  when  the  corresponding  search  region  is  not  visited,  which 
makes  it  possible  to  remove  the  nonlinear  interactions  in  the  objective  function  (3a),  yielding  (11a). 
It  is  well  known  that  unnecessarily  large  Big-M  values  lead  to  poor  continuous  relaxations  and 
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ultimately  slow  down  computation  time  [5].  In  the  case  of  the  SSP,  since  target  motion  is  linear, 
we  can  compute  these  values  based  on  the  maximum  distance  between  each  pair  of  targets. 

3.2  Partial  Path  Relaxations 

For  any  two  target  search  regions  i  and  j  (home  base  possibly  being  one  of  them),  the  following 
convex  NLP  produces  the  minimum  travel  distance  between  them. 

Problem  D 


Stj  =  mtin  !  I  (pi  +  (t~  Ti)Ui)  ~  (pj  +  (f  -  Tj)Uj)  1 1 

s.t.  max{r™’1,  r"un}  <t<  min{r™ax,  T™ax} 

We  let  SjtH+ 1  =  P_1h*n+1  be  the  minimum  travel  time  resource  expenditure  between  search 
region  j  and  the  recovery  location,  as  described  in  Section  2.  We  proceed  under  the  assumption 
that  the  home  station  and  the  recovery  location  are  the  same  physical  location,  therefore  Soj  = 
Sj:H+ 1  =  V^Sqj  =  P-1(5*n+1.  This  is  usually  the  case  in  the  type  search  planning  problems  we 
consider  and  it  imposes  no  limitations  on  onr  model  or  solution  procedures.  The  SSP  instance  of 
RP(O)  is  obtained  when  no  path  x  is  specified.  We  force  a  to  take  on  lower  bound  values  in  order 
to  allow  d  to  take  on  highest  possible  values.  When  this  is  done  a  and  t  can  be  eliminated  from 
the  problem,  resulting  in  the  following  NLP  in  the  search  time  d. 

Problem  RSSP(O) 


max 

d 

S.t. 


S  Qj  (1  —  exp  {—ajdj}) 

j€N 

dj  <  T™ax  -  T™n ,  Vj  G  N 

E  dj  <  T  -  2  min{5j  n+i} 

jeN  ''  v 


E  dj  <  D  -  2  min{  T,  ,n+l} 

jeN 

dn+ 1  0 


dj  >0  Vj  G  N 


For  any  partial  path  xe,  we  have  the  following  relaxed  NLP  as  a  special  case  of  RP(£b). 
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Problem  RSSP(£b) 


max  Qj  (1  —  exP  {— 1 oijdj}) 

a.d.t  .  ~ 

jeN 

s.t.  (II Pi  +  (a*  +  di-  Ti)ui  -  Pj  -  (a*  +  dt  +  Uj  -  Tj)uj ||  -  VUj)  <  0,  V(i,  j)  G  A(xe) 

0 at  +  di  +  titj  -  a,-)  <  0,  V(i,  j)  G  A{xt) 

dj  +  tij  <  T  —  (1  —  i^)(5qn+i 

j£N  (i,j)eA(xt) 

'Edj  +  Y2  tij  <  D  —  (1  —  If)Sftn+\ 

j£N  (i,j)eA(xf.) 

hdj  =  0,  Vj  G  N\N(xt) 

Uj>  0,  V(i,j)eA(xe) 

(9f)  through  (9j) 

SSP  can  be  solved  by  Algorithm  B&B  using  RSSP(ib)  relaxations  and  the  lower  bound  ini¬ 
tialization  heuristic  described  next. 

3.3  SSP  Heuristic 

Denote  the  optimal  solution  to  SSP  as  Z*,  we  observe  that  if  e  =  0  and  the  initial  guess  x  is 
provided  to  Algorithm  B&B  where  Z(x)*  =  Z* ,  then  the  number  of  NLP  solutions  required  to 
prove  x  =  x*  is  constant  regardless  of  how  branching  is  done  in  step  2.  This  is  a  direct  consequence 
of  the  fact  that  fathoming  only  depends  on  the  lower  bound.  Of  course  this  observation  is  not 
unique  to  our  problem  setting.  In  fact,  it  is  true  of  any  branch-and-bound  algorithm  provided  the 
algorithm  does  not  include  more  sophisticated  fathoming  rules.  This  observation  is  the  impetus 
to  develop  a  reliable  way  of  providing  initial  solutions  to  Algorithm  B&B,  possibly  eliminating  the 
need  for  complex  branching  strategies. 

In  order  to  provide  a  good  initial  solution  to  Algorithm  B&B,  we  propose  a  five-phase  heuristic 
that  relics  on  the  knowledge  that  solving  a  TOP-RDR  entails  finding  an  acceptable  balance  between 
dwell  and  transit  resource  expenditure.  Ramesh  and  Brown  [18]  outline  a  four-phase  heuristic  for 
the  GOP  using  a  bang-for-buck  ratio  that  relates  the  reward  at  each  node  to  the  bounds  on  transit 
time.  We  use  a  similar  idea  here,  however  since  rewards  and  transit  times  are  generally  not  known 
quantities,  we  consider  a  bang-for-buck  ratio  that  relates  expected  search  value  to  the  area  of  the 
search  region.  We  also  add  considerations  for  transit  time  by  clustering  targets  based  on  temporal 
and  spatial  proximity.  Throughout,  we  use  SSP(cc)  to  quickly  determine  the  value  of  search  plan 
at  each  iteration  in  the  heuristic. 

Phase  I  of  the  SSP  heuristic  initializes  parameters  and  solves  RSSP(O)  in  order  to  rule  out 
targets  that  have  low  search  value.  In  Phase  II,  remaining  targets  are  grouped  into  spatial  clusters 
based  on  geographical  region,  then,  using  a  clustering  parameter  A,  further  grouped  into  clusters 
based  on  their  earliest  arrival  times  rj"aa:.  In  our  SSP  application,  the  searcher’s  home  station 
is  generally  closer  to  target  arrival  locations  than  departure  locations.  Therefore  search  is  more 
likely  to  take  place  at  the  end  of  the  target’s  movement  track  than  it  is  at  the  beginning.  Earlier 
target  arrivals  represent  search  opportunities  that  are  eliminated  earliest  in  the  scenario.  Thus, 
earlier  arriving  targets  would  likely  be  searched  first  if  they  are  searched  at  all.  Targets  are  then 
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ordered  within  each  cluster  and  entire  clusters  are  ordered  with  respected  to  one  another.  Phase 
III  performs  a  feasibility  check  and,  if  necessary,  systematically  removes  the  lowest  value  targets  (of 
those  remaining)  from  consideration  in  order  to  find  an  initial  feasible  solution.  In  Phase  IV,  pairs 
of  consecutively  ordered  targets  that  belong  to  different  clusters  are  examined  in  a  cluster  seam 
refinement  process.  Finally,  all  remaining  targets  are  considered  for  removal  in  Phase  V.  Further 
details  of  the  SSP  heuristic  are  provided  in  Appendix  A. 

In  the  worst  case  the  SSP  heuristic  requires  solving  2n+2  NLPs.  This  occurs  when  a  feasible  path 
is  found  in  Phase  II,  and  all  targets  are  considered  for  removal  in  Phases  IV  and  V.  In  this  situation 
each  target  occupies  its  own  cluster.  This  can  be  prevented  in  well  posed  problem  instances  where 
the  clustering  parameter  A  is  chosen  appropriately  with  respect  to  the  arrival  times  7-™“.  Since 
the  NLP  subproblems  can  be  solved  quickly,  approximately  1/10  of  a  second  for  large  problems, 
the  heuristic  is  quite  fast  even  in  the  worst  case. 

4  Numerical  Results 

We  consider  SSP  where  smugglers  move  in  a  northwesterly  direction  as  they  attempt  to  transport 
illicit  material  from  the  south.  In  this  scenario,  we  assume  smuggler  movement  occurs  through 
corridors  defined  by  coastal  strips  of  likely  departure  and  arrival  locations  as  depicted  by  the  abstract 
map  shown  in  Figure  2.  Observe  that,  even  though  the  dotted  lines  intersect  on  the  spatial  map, 
target  search  regions  rarely  overlap  in  space  and  time. 

We  assume  searcher  and  target  performance  data  that  is  consistent  with  known  planning  factors 
for  P-3  aircraft  and  GO  FAST  smuggler  boats.  We  assume  a  24-hour  mission  execution  period 
(D).  Departure  time  uncertainty  data  is  randomly  generated  within  the  mission  planning  period 
with  uncertainty  ranging  up  to  four  hours.  Location  uncertainty  data  is  randomly  generated  where 
smugglers  are  equally  likely  to  depart  and  arrive  anywhere  on  the  aforementioned  coastal  strips. 
Expected  value  of  detecting  each  target  is  randomly  generated  within  the  range  [500,  5000]  lbs, 
corresponding  to  estimated  payload  capacity  of  GO  FAST  boats. 

We  solve  100  randomly  generated  problem  instances  with  3,  5,  7,  8,  and  10  smugglers  and 
compare  model  performance  using  Algorithm  B&B  with  heuristic  initialization  applied  to  SSP 
versus  solving  SSPM  directly  using  two  MINLP  solvers.  For  each  set  of  problem  instances  we 
deem  the  best  solver  to  be  the  one  that  identifies  the  optimal  solution  in  the  shortest  amount 
of  average  computing  time.  For  the  purposes  of  this  numerical  experiment,  Algorithm  B&B  is 
implemented  with  a  depth- first-search  strategy  and  the  optimality  tolerance  is  set  to  zero,  e  =  0. 
The  SSP  heuristic  is  implemented  with  temporal  clustering  parameter  A  =  6  hours  and  spatial 
clustering  corresponding  to  the  southwest  (bottom-left)  and  northeast  (top-right)  regions  shown 
in  Figure  2.  All  computations  are  done  on  a  64-bit  Windows  7  desktop  computer  (2x  Intel  Xeon 
3.46GHz;  RAM  24GB)  using  GAMS  23.8.  We  use  MINOS  to  solve  all  NLP  subproblems.  In  initial 
testing,  DICOPT  and  BONMIN  with  ECP  solver  option  appeared  to  be  the  most  effective  GAMS- 
based  solvers  for  directly  solving  SSPM.  Accordingly  we  limit  our  MINLP  numerical  results  to 
these  two  solvers.  For  brevity,  we  refer  to  Algorithm  B&B  with  heuristic  initialization  as  B&B 
and  BONMIN  with  ECP  solver  option  as  BONMIN(ECP).  For  the  remainder  of  this  paper,  unless 
otherwise  stated,  searcher  and  target  data  is  presented  in  nautical  miles,  nautical  miles  per  hour, 
hours,  and  pounds;  computation  runtimes  are  given  in  seconds;  and  optimality  gaps  are  reported 
as  a  percent  difference  from  the  optimal  objective  function  value. 
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Figure  2:  Target  movement  tracks  are  randomly  generated  with  origin  and  destination  points  chosen 
within  coastal  strips  marked  by  thick  solid  bars.  Given  an  origin  and  a  destination,  the  target  track 
goes  along  a  straight  line  within  corridors  indicated  by  dotted  lines  in  a  northwesterly  direction. 
The  searcher’s  home  station  is  identified  by  “x”.  Dotted  lines  shown  illustrate  possible  target 
movement  tracks.  Randomly  generated  target  movement  tracks  are  not  limited  to  those  depicted 
here,  but  stay  within  the  envelope  boundaries. 

As  indicated  in  Table  3,  all  three  solvers  are  able  to  solve  all  100  of  the  5-target  SSPs  to 
optimality  within  13  seconds  computing  time.  For  these  problem  instances  BONMIN(ECP)  is  the 
best  solver,  while  B&B  yields  the  slowest  average  runtime. 

We  use  performance  profiles  [9]  as  a  method  for  comparing  solver  runtimes.  Performance  profiles 
require  two  components:  performance  ratios  and  performance  metrics.  A  performance  ratio  is  a 
ratio  of  the  runtime  for  solver  s  on  problem  p  to  the  fastest  runtime  for  all  solvers  tested  on  problem 
p.  A  performance  metric  is  the  empirical  probability,  across  all  problems  p,  that  the  runtime  for 
solver  s  is  within  a  factor  of  k  of  the  fastest  solver  runtime.  A  performance  profile  is  a  distribution 
function  of  the  performance  metric  over  factors  k. 

We  see  in  Figure  3  that  BONMIN(ECP)  preforms  well  on  5-target  SSPs,  with  the  fastest  runtime 
for  nearly  90%  of  these  problems.  DICOPT  runtimes  stay  within  a  factor  of  3  of  the  fastest  runtime 
for  90%  of  problems.  B&B  lags  behind  the  MINLP  solvers,  with  runtimes  within  a  factor  of  7  of  the 
fastest  runtimes  for  approximately  80%  of  problems.  All  of  the  problem  instances  being  examined 
here  are  solved  in  13  seconds  or  less.  On  a  relative  (performance  profile)  scale  BONMIN(ECP) 
seems  to  be  the  clear  winner,  but  all  of  these  solvers  would  be  acceptable  to  planners  in  a  practical 
sense. 


16 


B&B 

BONMIN(ECP) 

DICOPT 

Num  Solved 

100 

100 

100 

Runtime  (sec) 

Average 

5.77 

1.25 

2.03 

St  Dev 

2.52 

0.30 

0.89 

St  Error 

0.25 

0.03 

0.09 

Median 

5.44 

1.22 

1.92 

90th  Percentile 

9.21 

1.70 

3.08 

Min 

1.54 

0.61 

0.61 

Max 

13.00 

1.98 

4.99 

Table  3:  Runtime  summary  for  5-target  SSPs.  Num  Solved  refers  to  the  number  of  problems  out 
of  100  that  were  solved  to  optimality  within  30  minutes.  BONMIN(ECP)  dominates  in  all  metric 
categories  and  B&B  yields  the  slowest  average  runtime. 


In  the  7  to  8-target  range  the  relative  performance  of  these  solvers  changes  dramatically.  Table  4 
and  the  performance  profile  plot,  Figure  4,  highlight  that  runtimes  of  BONMIN(ECP)  and  B&B 
on  7-target  SSP  instances  are  nearly  identical.  DICOPT  yields  the  slowest  runtimes  of  the  three 
solvers  tested. 


B&B 

BONMIN(ECP) 

DICOPT 

Num  Solved 

100 

100 

100 

Runtime  (sec) 

Average 

32.84 

30.47 

115.39 

St  Dev 

20.48 

16.31 

97.76 

St  Error 

2.05 

1.63 

9.78 

Median 

27.82 

26.19 

86.20 

90th  Percentile 

53.06 

52.99 

240.25 

Min 

8.77 

8.02 

15.90 

Max 

141.07 

81.96 

618.70 

Table  4:  Runtime  summary  for  7-target  SSPs.  BONMIN(ECP)  and  B&B  have  nearly  identical 
runtimes.  DICOPT  yields  the  slowest  runtimes  of  the  three  solvers  tested. 

Table  5  shows  that  B&B  is  the  best  solver  for  the  larger  8-target  problem  instances.  B&B 
runtimes  are,  on  average,  90  seconds  (1.5  minutes),  while  BONMIN(ECP)  runtimes  are  much  larger 
at  267  seconds  (4.5  minutes).  Observe  that  since  the  limiting  distribution  of  the  sample  mean  is 
normal,  and  considering  the  standard  error  of  the  average  runtimes,  we  can  say  with  high  (>  99%) 
confidence  that  the  true  mean  runtimes  for  B&B  on  all  8-target  problems  in  this  sample  space  are 
faster  than  that  of  the  other  two  solvers. 

The  performance  profile  plot  (Figure  5)  demonstrates  that  runtimes  for  BONMIN(ECP)  and 
DICOPT  in  nearly  all  problem  instances  are  several  times  larger  than  that  of  B&B,  with  over  half 
of  the  probability  mass  for  BONMIN(ECP)  being  in  the  k  =  3  to  k  =  9  range. 

When  considering  larger,  10-target,  SSPs  it  is  clear  that  B&B  is  the  only  viable  algorithm  among 
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Performance  Profile  for  5-Target  SSPs 


Figure  3:  Performance  profile  for  5-target  SSPs.  BONMIN(ECP)  preforms  well  on  5-target  SSPs, 
with  the  fastest  runtime  for  nearly  90%  of  these  problems.  B&B  lags  behind  the  MINLP  solvers, 
with  runtimes  within  a  factor  of  7  of  the  fastest  runtimes  for  approximately  80%  of  problems. 

the  three  solvers  tested.  Table  6  highlights  that  B&B  is  able  to  solve  97  out  of  100  problem  instances 
within  30  minutes  of  computing  time.  The  average  runtime  is  8.5  minutes  and  90%  of  the  problems 
are  solved  within  17  minutes  of  computing  time.  Meanwhile  BONMIN(ECP)  and  DICOPT  are 
unable  to  solve  any  of  the  SSP  test  problems  within  30  minutes. 

We  now  quantify  how  far  the  BONMIN(ECP)  and  DICOPT  solutions  are  from  the  10-target  SSP 
optimal  solutions  by  looking  at  the  reported  optimality  gaps  upon  termination.  Table  7  highlights 
that  when  BONMIN(ECP)  and  DICOPT  terminate  upon  reaching  the  30  minute  time  limit,  the 
solution  available  is  usually  far  from  optimal.  On  average  solutions  are  off  by  a  factor  of  at  least 
3.5  (optimality  gaps  in  excess  of  350%).  For  over  half  of  the  problems  tested,  DICOPT  is  unable 
to  provide  a  bound  on  the  optimal  solution  because  the  initial  MIP  for  the  linearized  subproblem 
is  not  solved  within  30  minutes. 

The  trend  continues  for  larger  problems.  On  a  set  of  25  randomly  generated  15-target  problem 
instances  B&B  solves  each  problem  instance  to  optimality  in  4667.77  seconds  (1.30  hours)  on  aver¬ 
age,  solving  18  out  of  25  problem  instances  within  2  hours.  BONMIN(ECP)  is  unable  to  solve  any 
of  these  problem  instances  within  2  hours,  terminating  with  an  average  optimality  gap  of  949%. 

We  are  able  to  gain  some  insight  into  why  B&B  outperforms  the  MINLP  solvers  as  the  problem 
size  increases  by  examining  the  branch-and-bound  enumeration  tree  of  a  representative  problem 
instance.  We  consider  a  10-target  SSP  instance  that  is  solved  in  324  seconds,  near  the  median 
runtime.  In  order  to  isolate  the  efficiency  gained  by  the  heuristic,  we  solve  this  problem  with  no 
initial  solution  provided  as  well  as  with  heuristic  initialization.  We  note  that  a  10-target  SSP  results 
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Performance  Profile  for  7-Target  SSPs 
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Figure  4:  Performance  profile  for  7-target  SSPs.  BONMIN(ECP)  and  B&B  each  have  the  best 
runtimes  in  approximately  half  of  the  test  problems.  DICOPT  yields  the  slowest  runtimes  of  the 
three  solvers  tested. 

in  an  enumeration  tree  of  nearly  20  million  nodes,  spanning  12  levels  deep.  At  each  node  we  solve 
RSSP(ib).  Clearly,  Algorithm  B&B  visits  only  a  small  fraction  of  these  nodes  due  to  pruning.  Any 
path  through  the  tree  that  has  length  12  is  a  path  that  visits  all  target  search  regions.  We  use  the 
term  perceived  depth  to  refer  to  the  depth  of  visited  nodes  in  the  enumeration  tree.  If  the  average 
perceived  depth  of  the  tree  were  large,  it  would  be  tantamount  to  enumerating  all  possible  paths 
x  G  X.  Thus,  in  order  for  our  B&B  algorithm  to  perform  efficiently  we  need  that  the  perceived 
depth  of  the  tree  remain  relatively  small  for  large  problems.  This  is  possible  due  to  the  dwell- 
to-transit  resource  trade-off  that  takes  place  when  we  consider  extensions  to  partial  paths  in  the 
enumeration  tree.  Table  8  shows  the  number  of  nodes  at  each  level  of  the  tree  for  a  representative 
problem  instance  with  and  without  heuristic  initialization.  We  see  that  the  tree  is  explored  no  more 
than  9  levels  deep,  as  the  partial  path  relaxation  provided  by  RSSP(cc^)  encounters  the  optimal 
solution  bound  fairly  shallow  in  the  tree.  The  majority  of  the  nodes  are  visited  in  levels  5-7  of  the 
tree.  Considering  that  a  5-target  SSP,  solvable  in  only  a  few  seconds,  yields  an  enumeration  tree 
that  is  7  nodes  deep,  the  perceived  depth  of  the  tree  for  larger  10-target  SSP  is  shallow  relative  to 
its  problem  size. 

Table  8  also  shows  the  benefit  of  using  the  SSP  heuristic  to  determine  the  initial  guess  for  B&B. 
Runtime  increases  proportionally  to  the  number  of  required  NLP  solutions  in  the  enumeration  tree. 
Having  a  good  bound  on  the  optimal  solution  reduces  the  total  number  of  required  NLP  solutions 
by  a  factor  of  2.3. 

We  conclude  this  section  with  some  remarks  on  the  performance  of  the  heuristic  presented  in 
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B&B 

BONMIN(ECP) 

DICOPT 

Num  Solved 

100 

100 

82 

Runtime  (sec) 

Average 

90.01 

267.82 

876.74 

St  Dev 

72.13 

185.82 

578.37 

St  Error 

7.21 

18.58 

57.84 

Median 

67.38 

191.67 

737.59 

90th  Percentile 

181.18 

500.62 

1800 

Min 

8.74 

65.58 

57.69 

Max 

415.33 

946.74 

1800 

Table  5:  Runtime  summary  for  8-target  SSPs.  B&B  is  the  best  solver  for  these  problem  instances. 
B&B  runtimes  are,  on  average,  90  seconds  (1.5  minutes),  while  BONMIN(ECP)  runtimes  are  much 
larger  at  267  seconds  (4.5  minutes) 


B&B  BONMIN(ECP)  DICOPT 
Num  Solved  97  0  0 


Runtime  (sec) 


Average 

515.91 

- 

- 

St  Dev 

697.87 

- 

- 

St  Error 

69.79 

- 

- 

Median 

313.06 

- 

- 

90th  Percentile 

1006.96 

- 

- 

Min 

33.27 

- 

- 

Max 

5048.39 

- 

- 

Table  6:  Runtime  summary  for  10-target  SSPs.  B&B  is  the  only  viable  solver  for  these  problems. 
BONMIN(ECP)  and  DICOPT  are  unable  to  solve  any  of  these  problems  within  30  minutes  of 
computing  time. 

Section  3.3  with  respect  to  the  3,  5,  7,  8,  and  10-target  SSP  test  set.  The  heuristic  correctly  identifies 
223  optimal  solutions  out  of  500  total  SSPs  tested.  Table  9  shows  that  the  accuracy  of  the  heuristic 
tends  to  diminish  as  the  number  of  targets  increases.  However,  the  average  relative  optimality 
gap  remains  within  the  1-3%  range  throughout.  Therefore  while  the  accuracy  rate  in  finding  the 
optimal  solution  decreases,  the  heuristic  does  not  miss  by  too  wide  of  a  margin  on  average.  The 
heuristic  is  able  to  get  within  7%  of  optimality  in  at  least  90%  of  all  problems  tested.  In  all  cases, 
the  average  accuracy  rate  is  driven  by  one  to  three  poor  performing  problem  instances.  While  the 
average  accuracy  diminishes,  the  runtime  remains  fairly  constant.  It  is  at  or  below  half  a  second  for 
all  problems  tested.  This  is  consistent  with  the  worst  case  runtime  analysis  presented  in  Section  3.3. 
Observe  that,  comparing  Table  9  to  Table  6,  on  average  the  heuristic’s  90th  percentile  optimality 
gaps  for  10-target  problem  instances  are  smaller  than  the  optimality  gaps  for  BONMIN(ECP)  and 
DICOPT.  A  problem-by-problem  comparison  reveals  that  the  heuristic  solutions  have  a  smaller 
optimality  gap  in  all  problem  instances. 
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Performance  Profile  for  8-Target  SSPs 


l 


k 


Figure  5:  Performance  profile  for  8-target  SSPs.  B&B  yields  the  fastest  runtimes  for  nearly  all 
problem  instances.  BONMIN(ECP)  and  DICOPT  runtimes  are  at  least  three  times  larger  than 
that  of  B&B  for  over  60%  of  problems  tested. 

5  Smuggler  Search  Problem  Model  Enhancements 

SSP  serves  as  a  baseline  model  for  solving  many  interesting  problems  that  arise  in  search  planning. 
Real-world  search  planning  scenarios  typically  involve  developing  search  plans  for  multiple  searchers, 
while  accounting  for  complicated  smuggler  movement  tracks.  We  show  how  our  model  can  be 
enhanced  to  capture  these  difficult  planning  problems. 

5.1  Multiple  Searchers 

We  consider  a  search  planning  operation  where  s  €  S  searchers  are  available.  We  model  this 
planning  problem  as  a  GOP-RDR  on  a  searcher-expanded  network  Gs  =  (Ns,  kls),  where  the  nodes 
are  sear cher-tar get  pairs  Ns  =  {(s,  j)  :  s  e  S,j  G  N}  and  the  arcs  As  =  {(s,i,j)  :  s  G  S,  ( i,j )  G  A } 
represent  the  transit  of  searcher  s  between  search  region  i  and  search  region  j.  Utilizing  the  vector 
forms  of  aj,dj  and  in  P.  we  allow  each  of  these  resource  variables  to  have  \S\  elements.  We 
denote  by  asj  and  dsj  the  respective  arrival  time  and  dwell  time  of  searcher  s  in  search  region  j, 
and  we  denote  by  tSjj  the  transit  time  of  searcher  s  from  search  region  %  to  search  region  j.  The 
multiple  searcher  SSP  is  stated  as  follows. 
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BONMIN(ECP) 

DICOPT 

Num  Bounds 

100 

44 

Gap  (%) 

Average 

504 

358 

St  Dev 

323 

223 

Median 

433 

316 

10th  Percentile 

236 

126 

90th  Percentile 

903 

551 

Min 

74 

17 

Max 

1893 

1182 

Table  7:  Optimality  gap  summary  for  BONMIN(ECP)  and  DICOPT  on  10-target  SSPs.  Num 
Bounds  refers  to  the  number  of  problem  instances  out  of  100  for  which  the  respective  solver  provided 
a  bound  on  the  objective  function  value  within  30  minutes  of  computing  time.  Reported  solutions 
generally  differ  from  the  optimal  solutions  by  a  large  margin,  upwards  of  350%  for  both  solvers. 
DICOPT  is  unable  to  report  an  optimality  gap  within  30  minutes  for  56  problem  instances. 


Problem  SSP-MS: 

max  I]  qj  (  1  -  exp  -J2asjdsjysj  >  )  (15a) 

a,d,t,x,y  V  l  ses  J  ) 

S-t.  (||Ps,i  T  (®s,i  T  dsi  Ts,i)^D,i  Ps,j 


& s,i  d s,i  ^ s,i,j  TsJ )USij 


•  •  •  ^  0, 

V(s,i,  j)  G  As 

(15b) 

^s,i  ds,i  &S, j')  %s,i,j  —  0? 

V(s,i,  j)  G  As 

(15c) 

xi  ds,j  +  x  <  Tsi 

jeNs 

Vs  G  S 

(15d) 

X  dsj  +  XI  <  Ds, 

jeNs  ( i,j)eAs 

Vs  G  S 

( 15e) 

n  .  >  rr-min 

as,j  —  '  s,j  ? 

V(s,  j)  G  Ns 

(15f) 

Cl  •  -1 -  (j  •  <C  q~rnfLX 

V(s,  j)  G  Ns 

(15g) 

&S,  0  0? 

Vs  G  A 

(15h) 

ds,n+ 1  0, 

Vs  G  A 

(15i) 

dsj  ^  0, 

V(s,  j)  G  Ns 

(18j) 

ts,i,j  ^  0, 

V(s,  i,j)  G  As 

(15k) 

Vs  r  xS) 

Vs  G  A 

(151) 

xs  G  X, 

Vs  G  A 

(15m) 

Each  expression  in  SSP-MS  is  a  direct  extension  of  its  SSP  counterpart  where  asj,dsj  and 
tSjij  are  the  arrival  time,  dwell  time,  and  the  transit  time  of  searcher  s  contained  in  the  vectors 
cij,dj  and  respectively.  We  allow  resources  T  and  D  in  (15d)  and  (15e)  respectively  to  vary  by 
searcher.  This  is  a  useful  feature  that  allows  the  model  to  account  for  heterogeneous  searchers.  In 
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With  Heuristic 

Without  Heuristic 

Level 

NLP  Solutions 

NLP  Solutions 

1 

1 

1 

2 

10 

10 

3 

100 

100 

4 

657 

702 

5 

2008 

3024 

6 

2051 

5061 

7 

714 

3354 

8 

85 

720 

9 

0 

40 

Total 

5626 

13012 

Runtime  (sec) 

324.78 

676.69 

Table  8:  Number  of  NLP  solutions  on  each  level  of  the  B&B  tree  for  a  representative  10-target 
SSP  instance.  The  perceived  depth  of  the  enumeration  tree  is  shallow  relative  to  its  problem  size, 
highlighting  the  resource  trade-off  that  motivates  the  GOP-RDR  and  highlights  the  benefit  of  the 
Algorithm  B&B. 

the  objective  function  (15a),  each  exponential  term  associated  with  search  region  j  in  the  random 
search  model  computes  detection  probability  by  accumulating  total  search  effort  for  all  searchers. 

Algorithm  B&B  can  be  used  to  solve  the  multiple  searcher  SSP.  We  modify  the  notation  in 
Algorithm  B&B  by  requiring  that  the  nodes  of  the  enumeration  tree  be  viewed  as  (s,j)~ pairs, 
where  s  G  S,j  G  N.  We  also  vectorize  R  and  8i>n+ 1  in  RP(R)  to  account  for  path  completion 
with  respect  to  each  searcher,  and  modify  the  path  completion  criterion  in  step  3  to  require  that 
(s,j)  =  (s,n  +  l),Vs  G  S.  In  principle,  this  can  be  done  for  an  arbitrary  number  of  searchers, 
however  the  size  of  the  enumeration  tree  is  exponential  in  the  number  of  searchers  <Sj .  Fortunately, 
real-world  applications  we  consider  call  for  search  planning  with  a  very  limited  number  of  searchers. 
Thus,  for  operational  and  computational  reasons  we  limit  our  numerical  results  to  2-searcher  SSPs. 

In  order  to  provide  a  good  initial  guess  to  Algorithm  B&B,  accounting  for  2  searchers,  we 
perform  our  SSP  heuristic  sequentially  as  follows.  We  set  the  temporal  clustering  parameter  A  =  6 
hours.  We  initialize  the  path  for  searcher  2  to  consist  only  of  arc  (0,  n  +  1).  This  ensures  a  feasible, 
but  certainly  not  optimal,  path  for  searcher  2.  We  then  run  the  (single-searcher)  SSP  heuristic 
for  searcher  1  and  fix  the  resulting  path.  We  then  do  the  same  for  searcher  2.  We  improve  the 
search  plan  by  considering  the  removal  of  targets  from  searcher  l’s  path,  performing  Phase  V  of 
the  SSP  heuristic.  We  then  do  the  same  for  searcher  2.  Lastly,  we  attempt  to  improve  the  plan  by 
allowing  the  searchers  to  swap  their  entire  search  paths.  The  modified  heuristic  returns  the  best 
search  plan  encountered  after  the  aforementioned  steps  are  completed.  We  do  not  consider  pairwise 
target  swaps  between  searchers  nor  do  we  consider  a  parallel  implementation  of  the  SSP  heuristic. 
Our  aim  is  to  quickly  provide  a  good  initial  solution  to  Algorithm  B&B.  While  we  expect  that  a 
parallel  heuristic  with  pairwise  target  swapping  would  produce  higher  quality  solutions  on  average, 
it  is  our  view  that  the  marginal  improvement  in  solution  quality  would  not  be  worth  the  increased 
computational  cost  and  complexity. 

To  illustrate  the  merits  of  the  SSP-MS  model,  we  consider  a  2-searcher,  10-target  SSP  example. 
As  the  numerical  results  in  Section  4  indicate,  a  significant  amount  of  computing  time  is  usually 
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Number  of  Targets 


3 

5 

7 

8 

10 

Avg  Time  (sec) 

0.19 

0.24 

0.30 

0.31 

0.34 

Min  Time  (sec) 

0.14 

0.17 

0.22 

0.19 

0.22 

Max  Time  (sec) 

0.40 

0.42 

0.65 

0.55 

0.67 

Num  Optimal 

80 

42 

41 

31 

29 

Avg  Gap  (%) 

1 

3 

2 

3 

2 

90th  Percentile  Gap  (%) 

3 

7 

6 

6 

6 

Max  Gap  (%) 

29 

34 

23 

48 

30 

Table  9:  SSP  Heuristic  performance  results.  The  heuristic  presented  in  Section  3.3  produces  run¬ 
times  that  are  at  or  below  the  half-second  mark  for  all  3,  5,  7,  8,  and  10-target  SSP  instances. 
Accuracy  in  finding  the  optimal  solution  appears  to  diminish  as  problem  size  increases,  however  the 
average,  90th  percentile,  and  maximum  optimality  gaps  remain  fairly  stable. 

required  to  solve  a  10-target  SSP.  When  another  searcher  is  considered  in  an  already  difficult 
problem,  squaring  the  number  of  possible  search  plans,  we  should  expect  that  the  B&B  runtime 
would  increase  substantially.  We  consider  search  planning  to  take  place  on  a  24-hour  cycle.  Thus, 
runtimes  that  exceed  24  hours  are  clearly  unacceptable.  It  is  apparent  that  solving  the  MINLP 
directly  using  either  BONMIN(ECP)  or  DICOPT  would  not  produce  a  solution  in  an  acceptable 
amount  of  time.  For  this  reason,  we  proceed  using  B&B. 

The  spatial  distribution  of  the  targets  in  this  example  is  split  (see  Figure  6),  with  4  of  them  in 
the  southwest  region  and  6  of  them  in  the  northeast  region.  The  speeds,  travel  times,  and  expected 
detection  value  of  these  targets  are  varied  (see  Table  11).  Given  two  aerial  search  assets,  a  P-3 
and  an  HC-130  (see  Table  10),  operating  from  the  same  home  station,  we  wish  to  develop  a  plan 
that  maximizes  the  expected  detection  value.  Complete  searcher  and  target  data  for  this  problem 
instance  are  provided  in  Table  10  and  Table  11  respectively. 


Searcher 

V 

V 

T 

D 

W 

A 

325 

205 

10 

24 

15 

B 

280 

240 

10 

24 

15 

Table  10:  2-searcher  SSP  example  data  is  based  on  known  performance  specifications  for  a  P-3 
(searcher  A)  and  a  HC-130  (searcher  B). 

We  select  this  particular  problem  out  of  the  100  randomly  generated  10-target  problem  instances 
because  it  exhibits  many  of  the  interesting  features  we  have  observed  in  our  study  of  these  problems. 
From  a  planner’s  perspective,  it  is  not  clear  how  to  form  a  good  plan  in  this  scenario,  let  alone 
an  optimal  one.  Targets  1,2,3,  and  5  have  high  expected  detection  values  q.  Good  judgment 
would  dictate  that  these  targets  would  likely  be  included  in  the  optimal  plan,  but  it  is  difficult  to 
determine  the  search  order,  the  searcher-to-target  assignment,  and  whether  or  not  other  targets  will 
be  searched  at  all. 

With  two  search  assets,  it  seems  reasonable  that  a  search  planner  would,  based  on  spatial  regions, 
assign  targets  2  and  3  to  one  searcher,  and  targets  1  and  5  to  the  other  searcher.  As  discussed  in 
Section  3.3,  we  assume  ascending  order  of  expected  arrival  times.  Furthermore,  since  the  targets  2 
and  3  appear  to  be  farther  away  from  the  home  station,  it  makes  sense  to  assign  the  faster  P-3  to 
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Figure  6:  10-target,  2-searcher  example  SSP  scenario  map.  4  targets  moving  in  the  southwest  region 
(lower  left)  and  6  in  the  northeast  region  (upper  right).  Movement  tracks  are  labeled  corresponding 
to  target  numbers  shown  in  Table  11. 

those  targets.  This  search  plan  is  reflected  in  Table  12  as  the  Manual  Plan  for  2  searchers. 

If  the  search  planner  only  has  one  search  asset  available,  developing  a  good  plan  somehow  appears 
even  more  difficult.  While  planning  to  search  all  the  high  value  targets  uses  a  substantial  amount  of 
travel  time,  which  is  a  limited  resource,  it  may  be  too  difficult  for  the  planner  to  determine  which 
of  these  targets  to  rule  out.  Therefore,  the  planner  may  assign  all  of  them  to  the  searcher,  again 
in  order  ascending  arrival  time.  This  search  plan  is  reflected  in  Table  12  as  the  Manual  Plan  for  1 
searcher. 

For  convenience,  we  define  a  p%-solution  to  be  a  solution  that  has  an  optimality  gap  of  p  percent. 
Table  12  shows  a  comparison  between  results  for  the  manual  plans,  search  plans  derived  from  the  SSP 
heuristic,  and  search  plans  computed  using  Algorithm  B&B.  The  manual  plan  for  a  single  searcher 
is  clearly  an  inferior  course  of  action,  a  17%-solution.  More  time  is  spent  transiting  between  search 
regions  than  is  spent  actually  searching  for  targets.  From  a  resource  trade-off  perspective,  we  can 
certainly  do  better.  The  single  searcher  heuristic,  a  3%-solution,  produces  a  plan  that  is  close  to 
optimal.  Two  of  the  high  value  targets  selected  (2  and  5)  are  correct  and  in  the  correct  order,  but 
the  heuristic  fails  to  recognize  the  value  of  adding  target  1  to  the  plan.  The  optimal  single  searcher 
plan  adds  target  1  because  the  additional  transit  resource  expenditure  is  worth  the  incremental 
bang-for-buck  in  searching  this  additional  target. 

The  manual  2-searcher  plan  has  an  optimality  gap  of  3%,  much  better  than  its  single  searcher 
counterpart.  This  is  possible  because  this  problem  instance  happens  to  have  target  movement  tracks 
for  high  value  targets,  two  in  each  region,  that  coincide  nicely  with  having  two  search  assets.  We 
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Target 

U 

r 

t 

P 

P 

P 

<1 

1 

59.5 

3 

4 

(24.9,  3.1) 

(9.1,  14.3) 

64 

4282 

2 

62.3 

9 

2 

(25.7,  10.1) 

(12.7,  20.2) 

55 

4554 

3 

61.7 

3 

2 

(27.0,  11.2) 

(16.3,  15.2) 

85 

3236 

4 

57.0 

10 

3 

(21.2,  1.9) 

(8.0,  14.7) 

59 

553 

5 

56.7 

4 

3 

(21.5,  2.6) 

(11.5,  13.2) 

27 

4687 

6 

59.2 

11 

1 

(26.6,  10.8) 

(13.0,  20.9) 

91 

937 

7 

55.3 

12 

1 

(21.1,  1.8) 

(5.8,  15.7) 

93 

661 

8 

64.5 

3 

4 

(27.7,  11.7) 

(12.9,  20.6) 

61 

837 

9 

64.4 

5 

4 

(24.4,  9.1) 

(14.3,  15.9) 

45 

1317 

10 

58.8 

7 

1 

(24.9,  9.5) 

(15.8,  15.4) 

66 

1290 

Table  11:  10-target  SSP  example  data  is  randomly  generated  as  described  in  Section  4. 


see  that  the  2-searcher  heuristic  obtains  its  plan,  as  expected,  by  starting  with  the  single  searcher 
heuristic  solution  and  adding  in  a  plan  for  the  second  searcher.  The  result  has  each  searcher  crossing 
between  geographic  regions  in  their  respective  paths.  This  plan  is  slightly  better  than  the  manual 
plan  and  yields  an  optimality  gap  of  2%. 


Expected  Search  Search  Transit 

Detection  Runtime  Order  Time  (hrs)  Time  (hrs) 


Manual  Plan,  1  searcher 

6514  <  0.1  sec  3-5-1-2  4.4  5.6 

Heuristic  Plan,  1  searcher 

7662  0.4  sec  5-2  6.5  3.5 

Optimal  Plan,  1  searcher 

7889  57.6  secs  1-5-2  6.4  3.6 

Manual  Plan,  2  searchers 

„„„  nl  A:  3-2  5.8  4.2 

12058  <  0.1  sec  &  ^  2  9 

Heuristic  Plan,  2  searchers 

..  A:  5-2  6.5  3.5 

12106  1.1  secs  B;  3_1  6A  3.6 

Optimal  Plan,  2  searchers 

,,,  A:  3-10-2  5.9  4.1 

12373  4.7  hrs  &  ^  ?  2  2  8 

Table  12:  Comparison  of  manual,  heuristic,  and  optimal  search  plans  for  10-Target  SSP  example 
with  1  and  2  searchers.  Search  order  includes  home  station  (not  shown)  at  the  beginning  and  end 
of  all  paths.  Plans  for  2  searcher  instances  are  labeled  A  and  B,  and  correspond  to  searcher  data 
in  Table  10.  Search  time  and  transit  time  refer  to  total  time  resource  expenditure  searching  and  in 
transit. 

The  optimal  2-searcher  plan,  depicted  in  Figure  7,  achieves  a  more  favorable  dwell-to-transit 
resource  trade-off  with  a  total  of  6.9  hours  spent  in  transit.  This  is  possible  because  of  two  not- 
so-obvious  modifications  to  the  manual  search  plan.  First,  target  10  is  added  to  searcher  A’s  path. 
This  enables  the  searcher  to  collect  reward  for  searching  another  target  while  allowing  target  2  to 
get  closer  to  the  searcher’s  location.  Thus,  more  time  is  spent  searching  than  in  transit.  Second, 
the  order  of  targets  1  and  5  are  swapped  in  searcher  B’s  path.  This  goes  counter  to  earliest  arrival 
time  ordering.  This  plan  takes  advantage  of  the  facts  that  the  movement  tracks  become  quite  close 
(where  the  dotted  lines  cross)  and  that  target  5  is  essentially  coming  straight  at  the  home  base  of 
the  searcher.  Combined,  these  modifications  yield  reduced  transit  time  and  increased  search  time. 
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Observe  that  in  the  optimal  2-searcher  plan  there  are  no  targets  searched  by  both  searchers,  see 
Figure  7.  One  might  conclude  that  solving  this  problem  is  equivalent  to  solving  two  separate  single 
searcher  SSPs.  This  is  not  the  case.  The  optimal  solutions  shown  in  Table  12  show  that  solutions 
are  not  nested  with  respect  to  adding  searchers.  While  the  optimal  1-searcher  plan  searches  regions 
1,  5,  and  2,  these  regions  are  not  allocated  to  the  same  searcher  in  the  optimal  2-searcher  plan. 

While  the  runtime  required  to  obtain  the  optimal  solution  is  relatively  long,  a  breath-first- 
search  of  the  first  three  levels  of  the  enumeration  tree,  done  in  9.6  minutes,  proves  that  the  objective 
function  value  of  the  heuristic  is  guaranteed  to  be  within  20%  of  optimality.  In  this  problem  instance 
a  search  planner  under  time  pressure  could  choose  to  use  the  heuristic  search  plan,  actually  (but 
unknown  to  the  planner)  a  2%-solution,  and  be  assured  that  this  plan  is  within  20%  of  optimality. 


Figure  7:  10-target,  2-searcher  example  SSP  optimal  solution  map.  Labels  indicate  search  plan 
ordering.  For  example,  label  “A2(10)”  is  the  searcher  A  to  target  10  assignment,  where  this  search 
region  appears  2nd  in  searcher  A’s  path.  All  search  regions  are  searched  in  the  bottom-right  to  top- 
left  direction  due  to  the  direction  of  target  motion.  The  size  of  each  rectangular  block  corresponds 
to  the  total  area  of  respective  search  region  during  the  time  when  the  searcher  is  performing  search 
actions  in  that  region. 

5.2  Complex  Target  Motion 

Real-world  scenarios  where  SSPs  arise  can  require  the  use  of  models  that  are  more  complex  than 
those  discussed  thus  far.  For  example,  a  particular  target’s  movement  track  may  not  be  along  a 
straight  line.  The  target  may  be  traveling  along  a  track  that  follows  a  particular  stretch  of  coastline, 
or  the  target  may  navigate  around  islands  or  other  geographic  obstacles.  It  is  also  possible  that  the 
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search  region  associated  with  a  target  changes  as  the  target  moves,  perhaps  due  to  changing  weather 
or  intelligence-driven  changes  to  the  uncertainty  ranges  themselves.  The  speed  of  the  target  may 
also  change  with  ocean  state  conditions  as  a  smuggler  travels.  All  of  these  considerations  can  be 
modeled  with  piecewise  linear  target  movement  tracks. 

We  consider  the  situation  where  target  motion  is  nonlinear,  but  can  be  approximated  by  piece- 
wise  linear  segments.  We  model  the  nodes  N  as  the  segmented  search  regions.  The  nodes  represent 
search  regions  as  in  SSP.  but  here  they  do  not  necessary  correspond  to  unique  targets.  Let  F  C  N 
be  the  set  of  first  segment  target  paths,  one  for  each  actual  target.  Let  B{i)  C  Nfi  e  F  denote  the 
set  of  search  region  segments  for  each  target.  The  piecewise  linear  target  motion  model  is  shown 
below. 

Problem  SSP-PWL: 


max  1  ~  exp  <  -  ajdjVj 

a,d,t,x,y  igF  y  [  i6SW 

s.t.  constraints  (9b)  through  (9m) 

SSP-PWL  only  differs  from  SSP  in  the  objective  function.  Similarly  to  SSP-MS,  we  sum  total 
search  effort  in  the  exponential.  In  this  case  the  summation  is  over  all  segments  in  the  piecewise 
linear  target  movement  track.  We  set  r?mm  and  r”ia:r  in  (9f)  and  (9g)  respectively  to  define  the 
connections  between  target  path  line  segments. 

Runtime  can  be  improved  by  reducing  the  size  of  the  arc  set  A,  thus  reducing  the  number  of 
partial  path  extensions  that  must  be  considered  in  step  2  of  Algorithm  B&B.  Clearly,  any  arc  (i,j) 
where  j  precedes  %  G  B(j)  should  be  eliminated  from  A.  Performing  this  elimination  procedure  for 
all  target  path  line  segments  reduces  the  dimension  of  X  in  Algorithm  B&B.  We  have  observed 
that  this  can  have  a  noticeable  effect  on  computation  time. 

To  illustrate  the  merits  of  the  SSP-PWL  model,  we  consider  a  2-searcher,  5-target  SSP  with 
piecewise  linear  target  tracks.  The  additional  segments  make  this  problem  equivalent  in  size  to 
a  12-target  SSP.  In  this  scenario,  one  of  the  targets,  target  1,  navigates  around  an  island  in  the 
southwest  region.  Once  past  the  island,  the  track  of  the  target  becomes  more  and  more  uncertain. 
Intelligence  estimates  are  solid  in  estimating  the  first  part  of  this  target’s  track,  but  analysts  are 
less  certain  about  the  target’s  arrival  location.  This  target  track  is  modeled  by  four  segments  (la, 
lb,  lc,  and  Id)  with  increasing  ranges  of  location  uncertainty.  Another  target,  target  2,  moves 
along  a  track  that  follows  the  coastline  in  the  northeast  region.  The  piecewise  linear  track  for  this 
target  is  captured  by  five  segments  (2a,  2b,  2c,  2d,  and  2e).  The  other  three  targets  move  along 
linear  tracks.  The  target  movement  tracks  for  this  scenario  are  shown  in  Figure  8. 

Target  data  for  this  scenario  is  given  in  Table  13.  Observe  that  the  location  uncertainty  p 
for  target  1  increases  by  25nm  in  each  subsequent  segment.  This  has  the  effect  of  increasing  the 
area  of  the  search  region  by  50  nm2.  The  expected  departure  times  r  in  subsequent  segments, 
j  e  N\F ,  are  derived  by  simply  propagating  the  first  segment  forward  in  time,  setting  r  to  the 
expected  arrival  time  of  the  previous  segment.  Target  departure  and  arrival  locations  are  chosen  to 
fit  the  scenario  description  and  to  make  the  map  illustration  readily  viewable.  All  other  values  are 
randomly  generated  as  described  in  Section  4. 

This  scenario  presents  a  difficult  challenge  for  a  search  planner.  With  2  search  assets  and 
5  targets  there  are  enough  resources  to  make  an  effort  to  search  all  targets,  but  it  is  not  clear 


Figure  8:  5-target  piecewise  linear  example  map.  Target  1  navigates  around  an  island  and  has 
increasing  range  of  location  uncertainty.  Target  2  moves  along  a  track  that  follows  the  coastline. 

how  to  order  the  search  or  how  determine  the  searcher-to-target  assignment.  Perhaps  the  most 
difficult  aspect  is  that  the  search  region  for  target  1  increases  in  size  as  it  gets  closer  to  the  home 
station.  Nominally,  a  good  search  plan  would  have  a  searcher  wait  at  home  station  while  the  target 
approaches,  but  doing  so  in  this  case  makes  it  harder  to  detect  the  target  later  in  the  mission 
execution  period  as  the  search  region  area  increases. 

The  optimal  solution,  shown  in  Figure  9,  is  obtained  in  10.5  hours  by  using  Algorithm  B&B  with 
settings  described  in  Section  5.1.  It  is  clear  from  the  map  that  the  path  for  searcher  A  achieves  a 
favorable  dwell-to-transit  resource  trade-off  as  search  regions  are  visited  at  a  point  that  seems  close 
to  the  home  station.  Searcher  B  is  forced  to  travel  a  relatively  long  distance  in  order  to  search  for 
target  1  at  a  point  where  its  search  region  has  a  small  area.  This  allows  target  5  to  be  searched  at 
a  point  that  is  close  to  the  home  station.  By  inspection,  this  optimal  search  plan  makes  sense,  but 
developing  it  without  using  optimization  techniques  entails  making  nonintuitive  assignments  in  the 
search  plan. 

A  breath- first-search  of  the  first  three  levels  of  the  enumeration  tree,  accomplished  in  2.7  hours, 
provides  a  solution  that  is  guaranteed  to  be  within  51%  of  optimality.  The  solution  turns  out  to 
be  a  1%-solution.  Given  only  a  couple  of  hours  for  planning,  a  search  planner  could  choose  to  use 
this  search  plan  in  lieu  of  the  heuristic  plan;  obtained  in  0.7  seconds,  but  with  no  optimality  gap 
guarantee.  The  heuristic  gives  a  53%-solution. 
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Target 

IJ 

T 

t 

P 

P 

P 

<1 

la 

64.1 

0 

2 

(23,  4) 

(10,  0) 

25 

2670 

lb 

64.1 

12.7 

2 

(10,  0) 

(10,  5) 

50 

2670 

lc 

64.1 

17.7 

2 

(10,  5) 

(10,  10) 

75 

2670 

Id 

64.1 

22.7 

2 

(10,  10) 

(10,  14) 

100 

2670 

2a 

61.4 

2 

3 

(28,  12) 

(21,  10) 

45 

890 

2b 

61.4 

9.1 

3 

(21,  10) 

(18,  10) 

45 

890 

2c 

61.4 

12.03 

3 

(18,  10) 

(17,  16) 

45 

890 

2d 

61.4 

18.0 

3 

(17,  16) 

(13,  17) 

45 

890 

2e 

61.4 

22.0 

3 

(13,  17) 

(13,  20) 

45 

890 

3 

62.7 

3 

1 

(26,  11) 

(15,  16) 

80 

2739 

4 

62.6 

10 

3 

(21,  1) 

(12,  13) 

62 

3914 

5 

55.5 

5 

1 

(28,  12) 

(13,  20) 

91 

2547 

Table  13:  5-target  piecewise  linear  SSP  example  data. 


Figure  9:  2-searcher,  5-target  piecewise  linear  optimal  solution  map.  Searcher  A  is  assigned  the 
path  3-5-2d,  searching  segment  d  of  target  2  last.  Searcher  B  is  assigned  the  path  lb-lc-4,  searching 
segments  b  and  c  of  target  1  first.  The  size  of  each  rectangular  block  corresponds  to  the  total  area 
of  respective  search  region  during  the  time  when  the  searcher  is  performing  search  actions  in  that 
region. 

6  Conclusions 

This  article  introduces  a  routing  problem  GOP-RDR  and  presents  a  specialized  branch-and-bound 
algorithm  that  is  built  upon  partial  path  relaxations  which  exploit  resource  trade-offs  that  are 
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inherent  in  these  problems.  We  formulate  a  search  problem  SSP  as  an  important  special  case  of  a 
GOP-RDR  and  provide  an  efficient  heuristic  for  computing  high  quality  solutions  in  a  small  amount 
of  time. 

Numerical  testing  on  randomly  generated  SSP  instances  reveals  that  our  branch-and-bound 
algorithm  outperforms  standard  MINLP  solvers  for  problems  with  7  or  more  targets.  In  large 
problem  instances,  with  10  targets,  our  algorithm  is  currently  the  only  viable  approach  to  obtain 
optimal  solutions  within  30  minutes  of  computing  time. 

We  propose  extensions  to  the  SSP  which  allow  practitioners  to  model  realistic  search  planning 
scenarios  that  involve  multiple  heterogeneous  searchers  and  complex  target  motion.  We  observe 
that  optimal  search  plans  can  usually  be  explained  by  an  intuitive  story,  but  obtaining  these  search 
plans  without  using  sophisticated  optimization  techniques  would  be  a  difficult  task. 
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A  SSP  Heuristic  Algorithm 

The  SSP  heuristic  begins  by  defining  A  as  the  temporal  clustering  parameter,  which  controls  how 
close  we  allow  target  clusters  to  be  with  respect  to  time.  We  assume  that  the  problem  of  interest 
can  be  separated  into  spatial  clusters  a  £  £  based  on  geographical  boundaries.  This  is  the  case 
in  our  SSP  example  where  smugglers  are  transiting  through  water  either  side  of  a  large  land  mass. 
Furthermore,  since  we  are  concerned  with  seagoing  smugglers  they  cannot  move  from  one  region  to 
another.  We  denote  by  Ka  the  set  of  targets  that  belong  to  spatial  cluster  a. 

SSP  Heuristic: 

Phase  I.  Initialization 

1.  Initialize  cluster  count  k  =  1,  order  index  o  —  1,  and  null  path  x  =  0.  Solve  RSSP(O)  and 
record  the  optimal  solutions  d*.  Construct  the  set  of  searchable  targets  N  =  N\{j  £  N  : 
d*  =  0}. 

2.  Compute  bang-for-buck  ratios  /3j  =  q3/ (fjpjUj),\/j  £  N. 

Phase  II.  Target  Clustering 

3.  For  each  spatial  cluster  a  £  £: 

Initialize  time  window  parameter  fa  =  min  While  fa  +  A  <  max  {r!nax}: 

jeNC\Ka  J  jGNCKv 

For  each  target  j  £  Ka: 

Assign  target  j  to  cluster  nk  if  r?maa:  £  [fa,  fa  +  A). 

Increment  fa  =  fa  +  A. 

If  one  or  more  targets  are  assigned  to  cluster  Kk  in  this  time  interval,  increment 
k  —  k  +  1. 
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4.  For  each  cluster  k 


Order  targets  j  G  Kk  in  ascending  value  of  Tjiax.  Compute  cluster  order  value  Uk  = 
niin{r"!a‘T  :  j  G  Kk}- 

5.  For  each  cluster  Kk,  considered  in  ascending  order  value  zyy 

For  each  target  j  G  Kk- 

Assign  search  order  Oj  =  o  and  increment  o  —  o  +  1. 

6.  Assign  order  O0  =  0  to  the  home  station  and  order  On+ 1  =  \N\  +  1  to  the  recovery  location. 

7.  Form  the  initial  path  x  by  setting  Xij  =  1  for  all  i  and  j  with  consecutive  orderings  0*  and 
Oj.  Solve  SSP(cc).  Save  incumbent  path  x  =  x.  If  the  problem  is  feasible,  save  the  heuristic 
objective  value  Z*H  as  the  optimal  objective  function  value  of  this  problem.  Otherwise,  set 
Z*H  =  —  oo. 

Phase  III.  Feasibility  Check 

8.  If  Z*H  =  — oo: 

For  all  targets  j  G  N,  considered  in  ascending  order  /3f. 

Do  Remove  j  procedure:  {  Remove  target  j  from  the  path  x  by  setting  xh]  =  0  (for 
i  :  Oi  =  Oj  —  1),  Xjji  =  0  (for  i!  :  Oy  =  Oj  +  1),  and  =  1  (for  i  :  Ot  —  03  —  1 
and  i!  :  O v  =  Oj  +  1).  Remove  j  from  the  set  N .  Solve  SSP(a;).  } 

If  a  feasible  solution  is  found,  set  Z*H  to  the  objective  function  value  of  this  solution 
and  go  to  step  9. 

Phase  IV.  Cluster  Seam  Refinement 

9.  Save  incumbent  path  x  =  x.  For  each  seam  between  clusters  Kk~\  and  Kk,  where  k  >  1: 

Let  j  be  the  target  in  that  last  order  position  in  cluster  Kk- 1-  Let  i!  be  the  target  in  the 
first  position  in  cluster  Kk ■  If  r™11  >  r™ax  and  (3j  <  fiy,  do  Remove  j  procedure  defined 
in  step  8.  If  Z*H  is  improved,  save  incumbent  path  x  =  x.  Otherwise,  reset  incumbent 
path  x  =  x. 

Phase  V.  Greedy  Target  Removal 

10.  For  each  target  j  G  N,  considered  in  ascending  order  [fi: 

Do  Remove  j  procedure  defined  in  step  8.  Solve  SSP(:r).  If  Z*H  is  improved,  save 
incumbent  path  x  =  x.  Otherwise,  reset  incumbent  path  x  =  x. 

11.  Return  heuristic  path  x  and  solution  Zb 
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